function [IRF,IRF_opt] = VARir(VAR,nsteps,ident,impact)
% =======================================================================
% Compute the IRFs for a VAR model estimated with VARmodel function
% =======================================================================
% [IRF IRF_opt] = VARir(VAR,nsteps,ident,impact)
% -----------------------------------------------------------------------
% INPUT
%   VAR     : structure, result of VARmodel function
%   nsteps  : number of nsteps to compute the IRFs
%
% OPTIONAL INPUT
%   ident   : type of identification: 'oir' (default) or 'bq'
%   impact  : size of the shock. 0 for 1stdev (default), 1 for unity
%
% OUPUT
%   IRF(t,j,k) : matrix with 't' steps, the IRF of the 'j' variable for the 'k' shock
%   IRF_opt
%     - IRF_opt.nsteps: number of steps
%     - IRF_opt.ident : identification chosen
%     - IRF_opt.impact: unit or 1 st dev
%     - IRF_opt.invA  : structural matrix such that invA*epsilon = u
% =======================================================================
% Ambrogio Cesa Bianchi, May 2012
% ambrogio.cesabianchi@gmail.com



%% Check inputs
%===============================================
if ~exist('impact','var')
    impact = 0;
end

if ~exist('ident','var')
    ident = 'oir';
end

%% Retrieve parameters and preallocate variables
%===============================================
c_case   = VAR.c_case;
nvars    = VAR.neqs;
nlags    = VAR.nlag;
beta     = VAR.beta;
sigma    = VAR.sigma;

IRF      = zeros(nsteps,nvars,nvars);


%% Compute the impulse response
%==============================

for mm=1:nvars

    % Compute the companion matrix (check for constant, trend, and/or 
    % exogenous variables and remove them
    switch c_case
        case 0
            beta_t = beta';
            F = [beta_t(:,1:nvars*nlags)
                eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];
        case 1
            beta_t = beta';
            F = [beta_t(:,2:nvars*nlags+1)
                eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];
         case 2
            beta_t = beta';
            F = [beta_t(:,3:nvars*nlags+2)
                eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];
    end
    
    Fn = eye(size(F,1));
      
    % Initialize the impulse response vector
    imp_res = zeros(nvars, nsteps);
    
    % Create the impulse vector
    impulse = zeros(nvars,1); 
    
    % Compute inv(A) for the impulse of the IRFs
    if strcmp(ident,'oir')
        [out chol_flag] = chol(sigma);
        invA = out';
        if chol_flag~=0
            disp('-------------------------------------');
            disp('Error: VCV is not positive definite');
            break
        end
    elseif strcmp(ident,'bq')
        Finf_big = inv(eye(length(F))-F); % from the companion
        Finf     = Finf_big(1:nvars,1:nvars);
        D        = chol(Finf*sigma*Finf')'; %Ident: u2 has no effect on y1 in the long run (this is different notation from notes!)
        invA = Finf\D;
    end
    
    if impact==0
        impulse(mm,1) = 1;  % 1 stdev shock
    elseif impact==1
        impulse(mm,1) = 1/invA(mm,mm); % 1 unity shock
    else
        disp('Error: impact must be either 0 or 1');
        break
    end

    % First period impulse response (=impulse vector)
    imp_res(:,1) = invA*impulse;
    big_impulse  = [(invA*impulse)' zeros(1, nvars*(nlags-1))]';
    
    % Recursive computation of impulse response
    for kk = 2:nsteps
        Fn = Fn * F; % this is the multiplier F^n
        big_imp_res   = Fn * big_impulse;
        imp_res(:,kk) = big_imp_res(1:nvars);
    end
    IRF(:,:,mm) = imp_res';
end

IRF_opt.nsteps    = nsteps;
IRF_opt.ident     = ident;
IRF_opt.impact    = impact;
IRF_opt.invA      = invA;
if strcmp(ident,'oir')
    IRF_opt.chol_flag = chol_flag;
else
    IRF_opt.chol_flag = 0;
end